Heat fluctuations in Ising models coupled with two difl"erent heat baths 
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Monte Carlo simulations of Ising models coupled to heat baths at two different temperatures are 
used to study a fluctuation relation for the heat exchanged between the two thermostats in a time r. 
Different kinetics (single-spin-flip or spin-exchange Kawasaki dynamics), transition rates (Glauber 
or Metropolis), and couplings between the system and the thermostats have been considered. In 
every case the fluctuation relation is verifled in the large r limit, both in the disordered and in the 
low temperature phase. Finite-r corrections are shown to obey a scaling behavior. 
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In equilibrium statistical mechanics expressions for the 
probability of different microstates, such as the Gibbs 
weight in the canonical ensemble, are the starting point 
of a successful theory which allows the description of a 
broad class of systems. A key-point of this approach is 
its generality. Specific aspects, such as, for instance, the 
kinetic rules or the details of the interactions with the 
external reservoirs are irrelevant for the properties of the 
equilibrium state. 

In non-equilibrium systems general expressions for 
probability distributions are not available; however, the 
recent proposal [ij, |3, Isl of relations governing the fluc- 
tuations is of great interest. They were formalized, for 
a class of dynamical systems, as a theorem for the en- 
tropy production in stationary states [3|. Fluctuation 
Relations (FRs) have been established successively for 
a broad class of stochastic and deterministic systems 
0, S, i, S S, S, E [iH (see [13 for recent reviews). 
They are expected to be relevant in nano- and biological 
sciences [13| at scales where typical thermal fluctuations 
are of the same magnitude of the external drivings. Test- 
ing the generality of FRs and the mechanisms of their 
occurrence in experiments IJ, ll5J or numerical simula- 



tions, particularly for interacting systems, is therefore an 
important issue in basic statistical mechanics and appli- 
cations. 

In this work we consider the case of non equilibrium 
steady states of systems in contact with two different 
heat baths at temperatures T„ {n = 1,2). In this case the 
FR, also known as Gallavotti-Cohen relation [3,] , connects 
the probability V{Q^'^\t)) to exchange the heat Q'"^'^) 
with the n-th reservoir in a time interval r, to that of 
exchanging the opposite quantity — Q'-"''(t), according to 



where 



In ^ ^ " = 0(")fr)A/3(") 

A/3(i) = (^ - ^) and A/3(2) ^ -A/^d). 



relation is expected to hold in the large r limit; in par- 
ticular r must be much larger than the relaxation times 
in the system. An explicit derivation of ([T]) for stochastic 
systems can be found in [l^]. In specific systems, the 



(1) 



This 



validity of ([T]) was shown for a chain of oscillators [17 1 
coupled at the extremities to two thermostats while the 
case of a brownian particle in contact with two reser- 
voirs has been studied in fl8|. A relation similar to 
([!]) has been also proved for the heat exchanged between 
two systems initially prepared in equilibrium at different 
temperatures and later put in contact [\v\. 

The purpose of this Letter is to study the relation ([T|), 
and the pre-asymptotic corrections at finite r, in Ising 
models in contact with two reservoirs, as a paradigmatic 
example of statistical systems with phase transitions. 
This issue was considered analytically in a mean field 
approximation in [20| where the distributions 7'(Q'"-'(t)) 
have been explicitly computed in the large r limit. Here 
we study numerically the model with short range inter- 
actions. This allows us to analyze the generality of the 
FR ([T]) with respect to details of the kinetics and of the 
interactions with the reservoirs, and to study the effects 
of finite r. We also investigate the interplay between 
ergodicity breaking and the FRs. 

We consider a two-dimensional Ising model defined by 
the Hamiltonian H — —J'YliUi) '^i'^j, where Ci = ±1 is a 
spin variable on a site i of a rectangular lattice with N = 
AI X L sites, and the sum is over all pairs (ij) of nearest 
neighbors. A generic evolution of the system is given by 
the sequence of configurations {cri(i), . . . , o'Ar(i)} where 
CTi(i) is the value of the spin variable at time t. We will 
study the behavior of the system in the stationary state. 
In order to see the effects of different kinetics, Monte- 
carlo single-spin-flip and spin-exchange (Kawasaki) dy- 
namics have been considered, corresponding to systems 
with non-conserved or conserved magnetization, respec- 
tively. Metropolis transition rates have been used; for 



single-spin-flip dynamics we also used Glauber transi- 
tion rates. 

Regarding the interactions with the reservoirs we have 
considered two different implementations. In the first, 
the system is statically divided into two halves. The left 
part (the first M/2 vertical lines) of the system interacts 
with the heat bath at temperature Ti while the right part 
is in contact with the reservoir at T2 > Ti. We have used 
both open or periodic boundary conditions. In the sec- 
ond implementation (for single-spin-flip only) each spin 
(Tj, at a given time t, is put dynamically in contact with 
one or the other reservoir depending on the (time depen- 
dent) value of hi = (1/2)| J2(j)i ^j\i where the sum runs 
over the nearest neighbor spins aj of Gi. Notice that hi 
is one half of the (absolute value) of the local field. In 
two dimensions, with periodic boundary conditions, the 
possible values of hi are hi = 0, 1, 2. At each time, spins 
with hi — 1 are connected to the bath at T = Ti and 
those with hi = 2 with the reservoir at T = T2. Namely, 
when a particular spin Ui is updated, the temperature 
Ti or T2 is entered into the transition rate according to 
the value of hi. Loose spins with hi = can flip back 
and forth regardless of temperature because these moves 
do not change the energy of the system. Then, as in 
the usual Ising model, they are associated to a temper- 
ature independent transition rate (equal to 1/2 or 1 for 
Glauber or Metropolis transition rates). This model was 
introduced in [2l[ and further studied in [24I. It is char- 
acterized by a line of critical points in the plane Ti,T2, 
separating a ferromagnetic from a paramagnetic phase 
analogously to the equilibrium Ising model. 

Denoting with tj^ the times at which an elementary 
move is attempted by coupling the system to the n-th 
reservoir, the heat released by the bath in a time window 
t C [s,s + t] is defined as 
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In the stationary state, the properties of Q'"''(t) will be 
computed by collecting the statistics over different sub- 
trajectories obtained by dividing a long history of length 
tp into many (tp /t) time-windows of length r, starting 
from different s. 
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FIG. 1: Upper panel (Ti, T2 > Tc): Heat PDs for Q*^' (on the 
left) and Q'^' (on the right) for a system with Ti = 2.9, T2 = 
3, size 10x10 and tp = 6-10* MCS. In the inset cr<-^'p(Q'^H'r)) 
is plotted against (Q'^'(^) - (Q*^'(^)))M^'- Curves for dif- 
ferent r collapse on a Gaussian mastercurve. Lower panel 
(ri,T2 < Tc): Same kind of plot for Ti = 1,T2 = 1.3, size 
10 X 10 and tp = 10^ MCS. The skewness of the distribu- 
tions (see text) is equal to 2.35, 0.589, 0.11 for the cases 
r = 4000, 10000, 40000 respectively. 



Ti, T2 > Tc. We begin our analysis with the study of 
the relation ([T]) in the case with both temperatures well 
above the critical value Tc — 2.269 of the equilibrium 
Ising model. In the following we will measure times in 
montecarlo steps (MCS) (1 MCS=N elementary moves). 

The typical behavior of the heat probability distri- 
bution (PD) is reported in the upper panel of Fig[T] 
for the case with static coupling to the baths, single- 
spin-fiip with Metropolis transition rates, Ti — 2.9, 
T2 = 3, and a square geometry with L = M = 10 
(Much larger sizes are not suitable because trajectories 



with a heat of opposite sign with respect to the av- 
erage value would be too rare). As expected, Q^^^(r) 
(Q(2)(r)) is on average negative (positive) and the 
relation < Q'^^^(t) > + < Q''^\t) >= is verified. 
Regarding the shape of the PD, due to the central limit 
theorem, one expects a gaussian behavior for r greater 
than the (microscopic) relaxation time (in this case 
it is of few MCS (- 5)), namely P(Q(")(t)) = 



(2^)-V2(a(")) 
(Q(")(r)) 



exp[ 
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y^T. This form is found 
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FIG. 2; Same parameters of Fig. [T] (different sizes). Upper 
Panel {Ti,T2 > Tc): ^{t, L) is plotted against t/L for differ- 
ent L. Inset; log [PiQ'^'''' {t))/V{-Q'^^''{t))] is plotted against 



(1/Ti - l/r2)Q(')C^)- Lower panel {Ti,T2 < T, 
ted against t/L for different L. 



.(2) 



is plot- 



with good accuracy, as shown in the inset of Fig. [1] 
where data collapse of the curves with different r 
is observed by plotting {ut Y^'^'P{Q'''^\t)) against 

(Q«(r)-(Q(^Hr)))M'^ 

In order to study the FR ([T]), we plot the logarithm 
of the ratio P(Q(")(t))/^(-Q^"H^)) as a function of 
A/3("')q("')(t), see Fig. [2] (inset). For every value of r 
the data are well consistent with a linear relationship 
(however, for large values of the heat the statistics be- 
comes poor), in agreement with the gaussian form of the 
PDs. To verify the FR ^, the slopes of the plot 



i:»(")(T) = 



, P(g'"'(r)) 
"^7'(-Q(")(r)) 

Q(")(t)A/3(") 



(3) 



must tend to 1 when r ^ oo. We show in Fig[2] the 
behavior of the distance e^"-' = 1 — D'^^\t) from the 
asymptotic behavior, for the case with static coupling to 
the baths. This quantity indeed goes to zero for large r 
(the same is found for dynamic couplings), e'"' depends 



in general on the temperatures, the geometry of the sys- 
tem and on r. Its behavior can be estimated on the basis 
of the following argument. 

In systems as the ones considered in this Letter, where 
generalized detail balance [16[ holds, the ratio between 
the probability of a trajectory in configuration space, 
given a certain intial condition, and its time reversed 
reads 



V{traj) 
Vi-traj) 



Q'"'(t)A/3("'-#^ 
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where AE — Q'^-'(t) + Q''-^^{t) is the difference between 
the energies of the final and initial state, and n' ^ n. 
For systems with bounded energy, Eq. @ is the starting 
point [16[ for obtaining the FR ([T]) in the large-r limit. 
Since Q'^'^^t) increases with r while AE is limited, in 
fact, the latter can be asymptotically neglected and, after 
averaging over the trajectories, the FR ^ is recovered. 
Keeping t finite, instead, from Eqs. (|3l4p one has that for 
the considered trajectory the distance of the slope from 
the asymptotic value is 
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We now assume that the behavior of e^"^ can be in- 
ferred by replacing AE and Q^"''(t) with their average 
values whose behavior can be estimated by scaling ar- 
guments. Starting with the average of Q*-"-'(r), we ar- 
gue that this quantity is proportional to Nfi^xT, Nfiux 
being the number of couples of nearest neighbor spins 
interacting with baths at different temperatures. This 
is because between these spins a neat heat flux occurs. 
In the model with static coupling to the baths one has 
Nfiux oc L. With dynamic coupling, instead, since every 
spin in the system can feel one or the other temperature, 
one has Nfi^x oc N. Concerning the average of AE, this 
is an extensive quantity proportional to the number N 
of spins. We then arrive at 



£(") ~ (r„,A/3("))- 



Lt-^ 

^-1 



static coupling to baths 
dynamic coupling to baths. 

(6) 
This result is expected to apply for sufhciently large r 
when our scaling approach holds. 

We remind that finite-r corrections have been shown 
to be of order 1 /r for the classes of dynamical systems 
considered in [3,l23|- The same is found in 2J] for models 
based on a Langevin equation [2^, in cases correspond- 
ing to the experimental setup of a resistor in parallel with 
a capacitor [15|. Faster decays (~ 1/r^) have been pre- 
dicted for other topologies of circuits 2J] . On the other 
hand, FRs in transient regimes, which are not considered 
in this Communication, are exact at any r (e = 0,Vr) 

In our model with static couplings the data of Fig. [2] 
confirm the prediction of the argument above: curves 



with different L collapse when plotted against x = t/L 
and e'"^ oc x~^ for sufficiently large r. Similar behav- 
iors have been found by varying the geometry (we also 
considered rectangular lattices with L > M), transition 
rates, and dynamics. The scaling prediction ^ has been 
verified for the Ising model with Kawasaki dynamics and 
squared geometry with L — 10,20,40, and in the case 
of the system dynamically coupled to the heat baths for 
sizes between L = 6 and L = 60. 



Ti,T2<Tc 

Let us first recall the behavior of the Ising model in 
contact with a single bath at temperature T < Tc- When 
N — oo the system is confined into one of the two pure 
states which can be distinguished by the sign of the mag- 
netization m{T) — {^/N)J2i^i'^i- This state is charac- 
terized by a microscopic relaxation time Teq (T) which is 
related to the fast flip of correlated spins into thermal 
islands with the typical size of the coherence length. At 
finite N, instead, genuine ergodicity breaking does not 
occur. The system still remains trapped into the basin 
of attraction of the pure states but only for a finite er- 
godic time Terg {N, T) , which diverges when N -^ oo or 
T — > 0. Then, as compared to the case T > Tc, there is 
the additional timescale Terg{N,T), beside Teq(T), which 
can become macroscopic. This whole phenomenology is 
refiected by the behavior of the autocorrelation function 
C{t - t') = (X;, cr^{t')ai{t)). When, for large iV, the two 
timescales Teq{T),Terg{N,T) are well separated, it first 
decays from C(0) = 1 to a plateau C{t - t') = m{T)'^ 
on times t — t' ~ Teq{T), due to the fast decorrelation 
of spins in thermal islands. The later decay from the 
plateau to zero, observed on a much larger timescale 
t ~ t' ~ Terg{N,T), signals the recovery of ergodicity. 
Notice also that, from the behavior of C{t — t') both the 
characteristic timescales can be extracted. 

The same picture applies qualitatively to the case of 
two subsystems in contact with two thermal baths, where 
each systems is trapped in states with broken symmetry 
for a time Terg{N,Tn) {Tf,rg{N,Ti) > Terg(iV, T2) since 
Ti < T2), which can be evaluated from the autocorrela- 
tion functions C^'^^t - t') = (E. CT|"''(^)^f''(i')), where 

(n) 

a\ denote spins in contact with the bath at T = T„. 
Since the FR is expected for r larger than the typical 
timescales of the system, it is interesting to study the 
role of the additional timescales Terg{N,Tn) on the FR. 
By varying Ti, T2 and N appropriately one can realize the 
limit of large r in the two cases with i) r <C Terg{N, T2) 
or ii) r ^ Terg{N, Ti). In case i), in the observation time- 
window T, the system is practically confined into broken 
symmetry states while in case ii) ergodicity is restored. 
Not surprisingly, in the latter case, we have observed a 
behavior very similar to that with Ti, T2 > Tc- The PDs 
for the case i) are shown in the lower panel of Fig. [T] The 
distributions are more narrow and non Gaussian. We 



calculated the skewness of the distributions defined as 
< {Q{t)- < Q{t) >f>/< (Q(r)- < Q{t) >)2 >3/2. 
This quantity is zero for gaussian distributions while for 
the case of the lower panel of Fig.[T]we found values differ- 
ent from zero which have been reported in the caption of 
the figure. This data show that the PDs slowly approach 
a Gaussian form increasing r. Moreover, differently than 
in the high temperature case, an asymmetry between the 
distributions of Q^^\t) and Q^'^''{t) can be observed. 

Regarding the slopes D("^(r), they converge to 1 also 
in this case even if the times required to reach the asymp- 
totic behavior are much longer than in the high temper- 
ature case (but always smaller than Terg(iV, T2)). This 
suggests that the FR ([T]) holds even in states where ergod- 
icity is broken and that the presence of the macroscopic 
timescales Terg{N,Tn) does not affect the validity of the 
FR in this system. Regarding the scaling of e'"-', the 
data are much more noisy than in the case Ti,T2 > T^ 
particularly for e^^' . Despite this, the data presented in 
the lower panel of Fig[2] for e*^^^ are consistent with the 
scaling ([6]) suggesting that it is correct also in this situ- 
ation. 

In this Letter we have considered different realizations 
of Ising models coupled to two heat baths at different 
temperature. We studied the fluctuation behavior of the 
heat exchanged with the thermostats and found that the 
FR (1) is asymptotically verified in all the cases consid- 
ered. We also analyzed the effects of finite time correc- 
tions and their scaling behavior. The picture is quali- 
tatively similar in the high and low temperature phases, 
although very different timescales are required to observe 
the asymptotic FR. 

The authors are grateful to A. Pelizzola and L. Ron- 
doni for useful discussions. 
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